Development and Validation of an Explainable Machine Learning Model to Assess the Prevalence Probability of Gastrointestinal Heat Retention Syndrome in Children: Cross-Sectional Study

Introduction

Human energy metabolism follows an energy balance equation, where energy intake equals total energy expenditure [1]. Energy intake occurs through dietary consumption, while total energy expenditure represents the sum of the organism’s daily basal metabolism (or resting energy expenditure), the thermic effect of food, and heat production from physical activity [2]. Basal metabolism refers to the energy expended at rest, accounting for approximately 60% of total daily energy expenditure. The thermic effect of food denotes the postprandial increase in metabolic rate, encompassing energy consumed for food processing and storage, as well as the metabolic effects of nutrient intake, constituting about 10% of total energy expenditure. Heat production from physical activity comprises both exercise and nonexercise activity thermogenesis. Conscious activities such as physical exercise constitute exercise thermogenesis, accounting for 0%-10% of total energy expenditure. Unconscious activities, including fidgeting and posture maintenance, represent nonexercise activity thermogenesis, contributing approximately 20% of total energy expenditure [3,4]. Since the Industrial Revolution began approximately 250 years ago, significant changes have occurred in human dietary patterns and lifestyles [5]. Decreasing food costs, advancements in processing technology, and increased food availability have progressively elevated energy intake in human diets. Concurrently, sedentary lifestyles characterized by physical inactivity and sun avoidance have gradually reduced energy expenditure. This energy imbalance ultimately results in a widespread energy surplus across populations, leading to various health problems [6,7]. Currently, most scholars posit that when energy intake exceeds total energy expenditure, the surplus is converted into fat for storage, consequently leading to obesity [8,9]. Our research team observed clinically that when energy intake surpasses the combined demands of basal metabolism, the thermic effect of food, heat production from physical activity, and energy expenditure for fat conversion, excess energy may be consumed through mechanisms such as enhanced metabolism or triggered inflammatory responses, manifesting as specific digestive tract symptoms. Animal studies have demonstrated that high-calorie diets can aggravate lipopolysaccharide-induced pneumonia through acetate-mediated macrophage polarization via the histone deacetylase 9/histone deacetylase 10-hypoxia-inducible factor-1α glycolysis axis. They also worsen 2,4-dinitrochlorobenzene–induced eczema through T helper 1/2 and T helper 17/regulatory T-cell imbalances, increased inflammatory factor expression, heightened oxidative stress, and alterations in the gut microbiota. Additionally, such diets elevate murine serum levels of the proinflammatory cytokines tumor necrosis factor-α and interleukin (IL)-6 while reducing anti-inflammatory IL-10. They upregulate Toll-like receptor 4, tumor necrosis factor-α, and IL-6 transcripts in the small intestine, alter intestinal permeability, and modify rat gut microbiota structure, inducing significant inflammatory changes in digestive tissues, thereby further substantiating this phenomenon [10-13].

Our team has provisionally termed this specific digestive tract syndrome—induced by enhanced metabolism and inflammatory responses resulting from high-calorie diets—as gastrointestinal heat retention syndrome (GHRS). GHRS primarily presents with gastrointestinal symptoms, including dry stool, straining during defecation, and strong fecal odor, and may be accompanied by dry mouth, increased eye discharge, epistaxis, oral ulcers, and elevated skin temperature of the palms and soles. Research demonstrates that GHRS is closely associated with childhood conditions such as recurrent respiratory infections, eczema, allergic purpura, and recurrent functional abdominal pain. Interventions targeting GHRS may therefore positively affect the prevention of such pediatric diseases [14].

The diagnosis of GHRS primarily relies on the diagnostic model and testing scale developed by our research team. The model is based on the Extreme Gradient Boosting (XGBoost) algorithm and achieved an internal validation accuracy of 93.03% [15]. The etiology of GHRS is closely associated with sociodemographic factors and lifestyle habits. Animal studies indicate that high-protein, high-calorie diets can induce GHRS. The established rat model exhibits characteristics including a preference for cooler areas, reduced bowel sounds, increased fecal characteristic scores, and elevated anal temperature. These changes are accompanied by elevated levels of the inflammatory factors IL-6, Iba-1 (ionized calcium-binding adapter molecule 1), and nuclear factor-kappa B p65, along with reduced levels of the blood-brain barrier structural proteins Claudin-5 and Occludin [16]. Our cross-sectional study in Beijing revealed variations in GHRS prevalence among children with different sociodemographic characteristics, including age, gender, household income, and parenting style. Children aged 3-5 years, girls, long-term urban residents, those with grandparents as primary caregivers, monthly household incomes of 10,000-20,000 CNY (1 CNY=US $0.15, as of June 19, 2026), antibiotic use before age 1 year, preterm birth, nonbreastfeeding, dinner after 7:00 PM, consumption of roasted, fried, or pan-fried foods 2 or more times per week, and frequent snack or dessert intake were identified as having a higher risk of GHRS. Conversely, daily indoor exercise exceeding 30 minutes and higher consumption of vegetables, fish, and soy products were identified as protective factors against GHRS [17]. Currently, research on GHRS is limited to the exploration of associated factors, and no prevalence prediction model has been developed; therefore, its practical clinical utility remains limited.

In recent years, machine learning (ML) techniques have been widely applied in disease screening. Zhu et al [18] used 8 demographic and clinical characteristics, including age, education level, hypertension, and a history of diabetes, to construct a prediction model for metabolic dysfunction–associated steatotic liver disease using 10 ML models, such as logistic regression (LR) and multilayer perceptron [18]. To enhance the credibility of model behavior and outcomes, ML models require explainability analyses. For instance, when Li et al [19] used the XGBoost model to predict prediabetes risk, they used the SHAP (Shapley Additive Explanations) method for explainability analysis. Their findings indicated that age, hypertension, and white blood cell count may be risk factors for prediabetes. When dealing with a large number of features, feature selection becomes necessary. Fitriyani et al [20] constructed an integrated prediction framework for cardiovascular disease by combining information gain feature selection, local outlier detection, and Bagging-Histogram Gradient Boosting. Through comparative testing against multiple traditional ML models, they demonstrated that feature selection combined with ensemble learning effectively enhances model discriminatory power and generalization capability [20].

This cross-sectional questionnaire survey investigated the positive rate of GHRS and its potential sociodemographic and lifestyle correlates among preschool children in Longgang District, Shenzhen. Using the collected data, we developed and validated an interpretable ML model to enable caregivers to estimate the probability that a child meets the study-specific GHRS scale criteria and subsequently developed a convenient assessment tool.

Methods

Data Collection

Data were collected through questionnaires from children attending 442 kindergartens in Longgang District, Shenzhen, between May and July 2021. The collected information encompassed basic demographic data, clinical symptoms and signs of GHRS, other symptoms and signs, and dietary behaviors. The inclusion criteria were (1) children enrolled in kindergartens in Longgang District during the survey period, and (2) guardians who were willing and able to complete the survey. The exclusion criterion was missing data on GHRS-related symptoms and signs, rendering confirmation of a GHRS diagnosis impossible.

Clinical Features

The variable data evaluated in this study were obtained through an online questionnaire survey completed by guardians based on the children’s actual circumstances. Children’s basic information encompassed 39 items, including kindergarten name, grade level, class, relationship to the informant, gender, date of birth, feeding method during the first 6 months of life, timing of animal protein complementary food introduction, and other relevant information. Clinical symptoms and accompanying signs of GHRS were 19 symptoms and signs observed within the preceding 2 weeks, such as increased palm and sole skin temperature, heat aversion, oral malodor, warm exhaled air from the mouth or nose, and thirst accompanied by a preference for cold drinks. Other symptoms and signs were derived from manifestations occurring within the preceding year and included listlessness, anorexia, sallow complexion, flaccid muscles, a preference for quietness, weak and feeble speech, reticence, excessive sweating, cold hands and feet, aversion to cold, dislike of raw and cold foods, and 51 additional symptoms or signs.

Research Results Evaluation

Symptoms and signs related to GHRS diagnosis were all included in the questionnaire. The questionnaire was completed by the children’s legal guardians, who possessed basic literacy skills. Following the questionnaire instructions, guardians provided responses based on their observations and judgments regarding the children’s symptoms and signs. This study used a guardian-administered questionnaire survey and did not implement blinding of investigators or outcome assessors. GHRS diagnosis was conducted using the “Self-Rating Diagnostic Scale for Gastrointestinal Heat Retention Syndrome in Children” developed by our team in prior research, which has demonstrated good reliability and validity [21,22]. The scale comprises 10 items with different weighting coefficients. Individuals with a cumulative score of 133 points or more were classified as GHRS positive, whereas those with scores below this threshold were classified as non-GHRS, as detailed in Table 1.

Data Preprocessing

Box plots were used to detect outliers in continuous variables. Outliers were defined as values above the upper quartile plus 1.5 times the IQR or below the lower quartile minus 1.5 times the IQR (Table S1 in Multimedia Appendix 1). Each outlier was replaced with the nearest boundary value to align with the main data distribution (Figure S1 in Multimedia Appendix 1). To maximize data validity, samples with missing values were excluded before model training. Categorical variables were predefined before model construction, and their data were processed using one-hot encoding. All sociodemographic subgroups were subjected to the same preprocessing procedures for outlier screening, missing value handling, and feature encoding, without differential treatment. Independent correlation analyses were performed to examine statistical associations between variables and GHRS and to assess the feasibility of model development. These analyses did not contribute to feature selection for the final model.

Variable Selection

In this study, GHRS diagnosis required the use of assessment scales. To prevent the trained ML models from learning the scale scoring rules directly or indirectly, we initially excluded variables related to GHRS clinical symptoms and accompanying signs (including scale items), retaining 467 other variables (Table S2 in Multimedia Appendix 1). Subsequently, univariate correlation screening and recursive feature elimination (RFE) were used for feature selection to enhance model performance and stability. Univariate correlation screening utilized the Pearson correlation coefficient to calculate the absolute correlation between each feature and the GHRS diagnostic outcome. Features with an absolute correlation coefficient of 0.1 or more were retained, whereas those exhibiting weaker correlations with the target variable were excluded to minimize interference from redundant variables. RFE is a widely used feature selection method in ML that eliminates unimportant features and ultimately identifies the optimal feature combination for achieving peak model performance. During the RFE process, 10-fold stratified cross-validation was applied, and features contributing minimally to model performance were recursively eliminated to derive the optimal feature subset for each model.

Sample Size Determination

According to the events per variable (EPV) rule, the minimum required sample size for the training set was calculated using the following formula:

Minimum required sample size for the training set = number of factors included × EPV/incidence of positive events in the training set

where the number of included factors refers to the number of features in the optimal feature subset obtained after univariate correlation screening and RFE. The positive event incidence rate in the training set represents the probability of a positive GHRS diagnosis. The EPV value was set at 10. According to existing research findings, the prevalence of GHRS among children is approximately 51.88% [17]. Therefore, the minimum required sample size for the training set should be at least 19.27 times the number of features in the optimal feature subset.

Model Development and Validation

Using stratified random sampling, the preprocessed dataset was divided into training and test sets at an 8:2 ratio. Final feature selection was conducted exclusively within the training set and cross-validation process to avoid information leakage from the test set. The test set was used solely for final model performance evaluation and did not participate in feature screening, resampling, hyperparameter tuning, or model selection. This study included samples only from multiple kindergartens in Longgang District, Shenzhen, and did not contain multiregional subsets; therefore, cross-subset heterogeneity was not assessed. After retaining complete cases, and considering that complete-case screening may alter the distribution of the training set, SMOTETomek (Synthetic Minority Oversampling Technique Tomek Links) was prespecified as a conservative resampling step in the training pipeline. Given the nearly balanced ratio of GHRS-positive to GHRS-negative cases, this step was not considered necessary to correct severe class imbalance but was retained as part of the primary analysis pipeline. SMOTETomek, feature standardization, and model training were integrated into an ImbPipeline (Imbalanced Learning Pipeline) and fitted exclusively on the training set and within cross-validation training folds. Validation folds and the test set were used solely for performance evaluation and were not involved in resampling, feature screening, hyperparameter tuning, or model selection, thereby minimizing the risk of information leakage. Seven ML models—Decision Tree, K-Nearest Neighbors, Linear Support Vector Machine (LinearSVC), LR, naive Bayes (NB), Random Forest (RF), and XGBoost—were used to estimate the probability of GHRS positivity. To optimize model performance, hyperparameter tuning was performed using grid search based on 3 rounds of 3-fold cross-validation to determine the final hyperparameters for each model using the optimal feature subset. Finally, because the outputs of the Decision Tree and LinearSVC models do not represent reliable probability estimates, Platt scaling (sigmoid) was applied for probability calibration to all models. Specifically, a calibration function was fitted using 3-fold cross-validation to convert model outputs into calibrated probability values ranging from 0 to 1.

After model construction, internal validation was performed. We conducted 10-fold stratified cross-validation on the training set. In each cross-validation iteration, feature standardization, SMOTETomek resampling, model training, and probability calibration were performed exclusively within the training folds, whereas the validation folds were used solely for performance evaluation. Performance metrics, including the area under the receiver operating characteristic curve (AUC), sensitivity, specificity, and Brier score, were calculated for each fold. The optimal threshold for each fold was determined using the Youden Index. Results were then aggregated across all folds to estimate the average performance on the training set and to derive the global optimal threshold. Using the optimal feature subset, optimal hyperparameters, and global optimal threshold obtained from the training set, the final model was trained. The model was subsequently validated on the test set, and performance metrics were calculated.

Model Performance Comparison

Model performance was compared across 3 dimensions: discriminatory power, calibration, and clinical utility. To evaluate discriminatory power, the AUC was primarily used to assess the model’s ability to distinguish between GHRS-positive and GHRS-negative samples. The AUC ranges from 0 to 1, with values closer to 1 indicating better discriminatory performance. Sensitivity and specificity were also calculated at both the fixed threshold (0.5) and the optimal threshold. Sensitivity reflects the model’s ability to identify positive samples, whereas specificity reflects its ability to identify negative samples. For calibration assessment, the Brier score was used to evaluate the accuracy of the model’s probability estimates. The Brier score ranges from 0 to 1, with values closer to 0 indicating better agreement between the predicted probabilities and the observed outcomes. Additionally, a calibration curve was plotted with the model-predicted probability on the x-axis and the observed positive rate on the y-axis to visually assess the model’s calibration performance. To further quantify calibration, the calibration slope and calibration intercept were calculated. The calibration slope typically approaches 1, with values closer to 1 indicating better agreement between the trend of the predicted probabilities and the observed positive rate. The calibration intercept typically approaches 0, with values closer to 0 indicating lower overall prediction bias. For clinical utility evaluation, decision curve analysis (DCA) was used to calculate the net benefit across different threshold probabilities. A higher net benefit indicates greater practical value of the model at a given clinical decision threshold. Decision curves were then plotted to visually compare the net benefits of each model with those of the “intervention for all” and “intervention for none” strategies. Finally, based on the performance of the above evaluation metrics in both the training and test sets, the optimal model was selected according to its overall practical utility.

Model Explanation

This study used the SHAP method to interpret the optimal ML model. SHAP is a game theory–based technique that ranks the importance of input features and explains model outputs, thereby addressing the limited interpretability of ML models. By calculating each feature’s contribution to the model output, SHAP provides both global and local interpretability, enhancing model transparency and explainability. In addition, a predicted probability distribution plot for the optimal model was generated to illustrate the agreement between model-predicted probabilities and actual outcomes, thereby supporting the reliability of the model evaluation results.

Network Calculator

To facilitate exploratory self-assessment at home, we integrated the optimized final model into a web application developed using the Streamlit framework. After users complete the online questionnaire and enter values for all predictive features, the application calculates and returns the probability that the child meets the study-defined GHRS scale criteria.

Statistics

The aforementioned ML models were developed using the scikit-learn and Keras libraries in Python (version 3.11; Python Foundation). Continuous variables with normal distributions are presented as mean (SD) and were compared using the unpaired (2-tailed) t test. Categorical variables were represented using one-hot encoding. Independent correlates in the overall cohort were identified through univariate LR analysis, collinearity testing, stepwise Wald regression analysis, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis, and multivariable LR analysis.

Ethical Considerations

This study protocol adhered to the principles of the Declaration of Helsinki and received approval from the Ethics Review Committee of Beijing University of Chinese Medicine, Shenzhen Hospital (Longgang; approval number SZLDH2021LSYM-03). As a cross-sectional study, this research was not registered in a clinical trial registry before data collection. This study was conducted with informed consent obtained from kindergarten administrators and the guardians of participating children. Data collection complied with the Ethical Review Measures for Life Science and Medical Research Involving Humans and the Personal Information Protection Law of the People’s Republic of China. The dataset was rigorously anonymized and deidentified to ensure that no participants could be identified, and no participant identities could be discerned from any images or multimedia appendices associated with this article. Participants received standardized compensation determined by the research team. Detailed compensation terms are not disclosed beyond the participants and members of the research team.

This study was reported in accordance with the TRIPOD+AI (Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis + Artificial Intelligence) guidelines. The complete TRIPOD+AI abstract and full reporting checklists are provided in Multimedia Appendices 2 and 3, respectively.

Results

Baseline Clinical Characteristics

Following questionnaire distribution, 120,198 questionnaires were returned. Of these, 11,751 were excluded because of missing data on GHRS symptoms and signs, leaving 108,447 questionnaires for analysis. The demographic characteristics of the children are presented in Table 2. As shown in Table 2, significant differences (P<.001) were observed between the GHRS-positive and GHRS-negative groups with respect to gender, grade level, age, height, and body weight, indicating that these demographic variables may contribute to estimating the probability that a child meets the GHRS scale criteria. Details of the study design are presented in Figure 1.

‎

Independent Correlative Factors

To determine the association between GHRS prevalence and questionnaire-collected variables while preliminarily assessing the feasibility of constructing a GHRS prevalence prediction model, we investigated independent correlates of GHRS in children using a dataset comprising 108,447 samples. After creating dummy variables and applying one-hot encoding to categorical variables, 469 variables were obtained. Variables such as “unclear,” which were designed to improve questionnaire completeness but lacked analytical significance, were excluded, leaving 467 variables. Univariate LR analysis identified 433 potential correlates associated with GHRS prevalence (P<.01; see Table S2 in Multimedia Appendix 1 for the complete list of P values). After collinearity testing, 335 potential correlates were retained (variance inflation factor [VIF]<10). Following stepwise regression with Wald tests, 170 potential correlates were identified (P<.001), and LASSO regression analysis identified 113 potential correlates (coefficient≠0). Finally, multivariable LR analysis identified 59 independent correlates associated with GHRS prevalence, including 10 protective factors and 49 risk factors (Table S2 in Multimedia Appendix 1 and Tables 3 and 4).

Data Preprocessing and Feature Selection

Using the dataset comprising 108,447 samples, records with missing values were excluded, yielding a complete-case sample of 49,798. Significant differences were observed in 304 out of 467 (65.1%) variables between the included and excluded pediatric groups (Table S3 in Multimedia Appendix 1). After stratified random sampling at an 8:2 ratio to create the training and test sets, the training set contained 39,838 samples, and the test set contained 9960 samples. Subsequently, the feature selection pipeline was reestablished. In the training set, after univariate correlation screening using the Pearson correlation coefficient, a total of 162 features were retained (Table S4 in Multimedia Appendix 1). Following RFE and a 10×10 cross-validation process, the optimal feature subset for each model was obtained (Table 5).

Sample Size Determination

Using the maximum number of optimal features among all models, the minimum required sample size for the training set was calculated to be 3122. With the training and test sets divided at an 8:2 ratio, the corresponding minimum required test set sample size was 625 cases. Therefore, the total minimum sample size required for model development and validation was 3747 cases.

Model Development and Validation

On the training set, we performed 10×10-fold internal cross-validation and hyperparameter grid search (Table S5 in Multimedia Appendix 1) to construct 7 ML models (Table 6, Figure 2). The results showed that XGBoost achieved the best AUC performance (mean 0.7417, SD 0.0086; 95% CI 0.7352-0.7482). Using 0.5 as the classification threshold, RF demonstrated the highest sensitivity (mean 0.6752, SD 0.0115; 95% CI 0.6666-0.6839), whereas NB exhibited the highest specificity (mean 0.7387, SD 0.0086; 95% CI 0.7323-0.7452). At the optimal threshold determined by the maximum Youden index, NB achieved the highest sensitivity (mean 0.6977, SD 0.0384; 95% CI 0.6688-0.7267), whereas LR demonstrated the highest specificity (mean 0.6783, SD 0.0304; 95% CI 0.6554-0.7013). For the Brier score, XGBoost exhibited the best performance (mean 0.206, SD 0.0031; 95% CI 0.2037-0.2084), with LR ranking second (mean 0.206, SD 0.0030; 95% CI 0.2037-0.2083).

On the test set, we evaluated the discriminatory performance, calibration, and clinical utility of each model (Table 7, Figure 3). Among the discrimination metrics, the highest AUC was achieved by XGBoost; the highest specificity at the 0.5 threshold by RF; the highest sensitivity at both the 0.5 threshold and the optimal threshold by LinearSVC; and the highest specificity at the optimal threshold by NB. DeLong test results indicated statistically significant differences in AUC among the models (Table 8). Among the calibration metrics, the mean Brier scores for XGBoost, LR, and RF were closest to 0.20, with RF exhibiting the narrowest CI for the Brier score. This finding suggests smaller fluctuations in calibration performance and slightly greater stability across different sampling scenarios. Regarding the calibration slope and calibration intercept, LR and XGBoost performed slightly better than RF, indicating that RF exhibited mildly conservative calibration and a tendency to underestimate risk relative to these models. For clinical utility assessment, bootstrapped CIs for the net benefit derived from DCA (Table S6 in Multimedia Appendix 1), together with the decision curves, showed that within the threshold probability range of 0.01-0.90, the net benefit of the RF model exceeded those of both the “treat-all” and “treat-none” strategies. Within the threshold probability range of 0.13-0.85, the net benefit of the RF model was substantially higher than that of the other models, demonstrating superior clinical decision-making utility across the study’s threshold spectrum.

The above results indicate that, among the 7 ML models, the RF model did not outperform all other models across every performance metric. XGBoost and LR demonstrated slight advantages in certain discrimination or calibration metrics. According to the prespecified comprehensive utility-oriented selection criterion, RF was chosen as the primary analysis model because of its narrower CI for the Brier score, acceptable calibration performance, and relatively greater net benefit in DCA, indicating superior overall utility.

‎

‎

Model Description

To improve the interpretability and credibility of the established prediction model, we used the SHAP method to analyze the contribution of each variable to the model output. We also generated a probability distribution plot for the RF model on the test set to further support the reliability of its prevalence probability predictions.

The SHAP method provides model interpretability through both global and local explanations. Global interpretability characterizes overall model behavior and is illustrated by SHAP summary bar plots, SHAP summary dot plots, and SHAP decision plots, whereas local interpretability is represented by SHAP dependence plots. Given the large number of features and their lengthy names in this study, we created a feature ID-feature name mapping table (Table S7 in Multimedia Appendix 1) to facilitate simplified feature representation. For global interpretation, SHAP summary bar plots were used to evaluate feature contributions to the model based on mean absolute SHAP values and to display the top 20 contributing features in descending order (Figure 4A). SHAP summary dot plots were used to visually demonstrate the direction and magnitude of each feature’s influence on model outputs (Figure 4B). SHAP decision plots illustrated the decision-making process for the first 100 samples (Figure 4C). The 5 most influential features were Feature19 (always quiet: none), Feature20 (speak in a low voice: none), Feature17 (soft muscles, not as strong as other children: rarely), Feature18 (soft muscles, not as strong as other children: sometimes), and Feature13 (yellow complexion, not as ruddy as other children: none). Higher values for Feature19 (always quiet: none) and Feature20 (speak in a low voice: none), together with lower values for Feature17 (soft muscles, not as strong as other children: rarely), Feature18 (soft muscles, not as strong as other children: sometimes), and Feature13 (yellow complexion, not as ruddy as other children: none), were associated with a higher predicted probability of GHRS in the model output. For local interpretability, SHAP dependence plots were generated for these 5 features to further illustrate their individual effects on model predictions (Figure 4D).

‎

The probability distribution plot was used to assess the agreement between the model-predicted sample distribution and the observed positive rate (Figure 5). Test-set samples were stratified into 3 subgroups according to their predicted probabilities of meeting the study-specific GHRS scale criteria: low probability (0%-30%), moderate probability (30%-60%), and high probability (60%-100%). The low-probability group (n=9960) consisted of 1898 (19.06%) test-set samples, including 311 positive cases, corresponding to an observed positive rate of 16.39%. This observed rate was substantially lower than the upper threshold of 30%, suggesting that the model was conservative within this probability range and could preliminarily identify children with a low likelihood of meeting the study-specific GHRS scale criteria. The moderate-probability group (n=9960) comprised 5302 (53.23%) test-set samples, including 2301 positive cases, corresponding to an observed positive rate of 43.41%. This finding suggests relatively good calibration performance within this probability range. The high-probability group (n=9960) included 2760 (27.71%) samples, of which 2008 were positive cases, corresponding to an observed positive rate of 72.75%. This result indicates that the RF model was able to appropriately stratify children with a high likelihood of meeting the study-specific GHRS scale criteria. Overall, the predicted probabilities of GHRS positivity generated by the RF model were consistent with the observed trend in positive rates, demonstrating satisfactory risk-stratification performance.

‎

Implementation of the Web Calculator

To facilitate public use, we integrated the RF model into a web application developed using the Streamlit framework for use in everyday settings. We established a definitive mapping between the 158 features included in the RF model and the questionnaire items used in the Streamlit application. The questionnaire content is consistent with the survey instrument used for data collection in this study, and all items can be understood and completed by a child’s guardian for the purpose of information collection (Table S8 in Multimedia Appendix 1). After the child’s guardian completes the electronic questionnaire consisting of 75 single-choice questions, the application automatically estimates the probability that the child currently meets the GHRS diagnostic criteria, outputs the corresponding disease probability category, and provides recommendations based on the results. This application can be accessed online [23].

Discussion

This study analyzed independent correlates of childhood GHRS using data from a cross-sectional survey, investigated and compared 7 ML models, and developed a prevalence prediction model for childhood GHRS using the model with the greatest overall utility. The model was subsequently deployed online, enabling its use as a tool for exploratory self-assessment of the probability that a child meets the GHRS scale criteria. This tool allows caregivers to estimate the likelihood of GHRS in children based on sociodemographic characteristics, lifestyle factors, and daily manifestations.

Analysis of independent correlates revealed that male gender; only-child status; a mother’s third pregnancy; regular lunch consumption; a maternal education level of junior high school, senior high school, or technical secondary school; absence of muscle softness; absence of excessive sweating during light activity; infrequent diarrhea after consumption of raw or cold foods; and lack of parental concern about picky eating were protective factors. According to formulas used to estimate energy requirements in children and adolescents, boys have significantly higher daily total energy expenditure than girls of the same body weight. Consequently, high-calorie diets may have a weaker impact on boys than on girls [24]. A cross-sectional study conducted in Shanghai reported that children with siblings experience more childhood trauma and adverse parenting experiences, whereas children without siblings receive more attentive care [25]. Third-born children, benefiting from their later birth order, often receive increased attention, leading to more refined dietary control. Children of mothers with secondary education demonstrated a lower likelihood of GHRS than those of highly educated mothers, contradicting findings from French and Dutch studies [26,27]. This divergence may stem from regional dietary and cultural differences. Parents’ lack of concern about their children’s picky eating and consistent lunch habits indicates appropriate dietary management. The absence of muscle flaccidity and excessive sweating during minimal activity suggests that the child’s energy expenditure falls within normal limits [28,29]. Infrequent diarrhea after consuming raw or cold foods indicates a healthy digestive state [30]. Conversely, a history of functional constipation, frequent or occasional thick tongue coating, frequent sallow complexion, frequent eye redness or susceptibility to gum problems, poor self-control with tantrums during discipline, consistently eating less than peers, persistent excitability, a history of functional dyspepsia, and low activity levels were identified as major risk factors. A history of functional constipation, consistently eating less than peers, frequent sallow complexion, and a history of functional dyspepsia indicate an unhealthy gastrointestinal state or increased susceptibility to adverse factors. Research indicates that a thick, greasy tongue coating is closely associated with metabolic disorders, and its microbiota may promote gastritis and influence the development of various chronic diseases [31]. Frequent eye redness and gum discomfort represent manifestations of the body’s inflammatory response, whereas mental agitation and poor self-control may reflect excessive energy expenditure associated with excessive energy intake, thereby serving as risk factors for GHRS [32]. Analysis of independent correlates revealed associations among sociodemographic factors, children’s daily behaviors, lifestyle status, and GHRS, supporting the feasibility of developing a prevalence probability assessment model using relevant variables. To protect respondents’ privacy and respect their preferences, all items in the online questionnaire were optional, resulting in a relatively high proportion of missing data. To assess whether directly removing samples with missing values would degrade model performance, this study constructed 5 imputed datasets using multiple imputation by chained equations and combined model evaluation metrics according to Rubin’s rules for sensitivity analysis. The results showed that, after multiple imputation, the RF model exhibited AUC, sensitivity, specificity, and Brier score values comparable to those of the original-case model on the test set. Although potential selection bias persisted, the impact of excluding samples with missing values on the model conclusions remained within an acceptable range (Table S9 in Multimedia Appendix 1). Subsequently, we compared baseline characteristics across all variables between children included in and excluded from the original-case analysis. The results indicated significant differences in 304 of the 467 (65.1%) variables between the included and excluded pediatric groups (P<.01; Table S3 in Multimedia Appendix 1). Furthermore, we noted that the vast majority of variables in this study were binary. Following dataset imputation, some imputed values produced logical inconsistencies. For example, both “father’s first child” and “father’s third child” could be imputed as “yes” for the same participant, rendering the postimputation data clinically implausible and potentially leading to misleading model conclusions. Therefore, to ensure data validity during preprocessing, this study excluded all samples with missing values. As complete-case analysis may alter the distribution of the training set, we included SMOTETomek as a conservative resampling step in the primary analysis pipeline. However, the ratio of GHRS-positive to GHRS-negative cases in this study was nearly balanced. Thus, SMOTETomek was not required to correct severe outcome-class imbalance, and synthetic sampling could potentially affect probability calibration. We therefore conducted an additional RF sensitivity analysis without SMOTETomek. The results showed minimal differences between the 2 pipelines in AUC, Brier score, calibration slope, calibration intercept, and DCA, indicating that the main analysis conclusions were robust to the inclusion or exclusion of this resampling step (Table S9 in Multimedia Appendix 1). As the weight assigned to a given variable during model training did not fully align with the results of the multivariable logistic analysis, and because the optimal feature set varied across model types, this study did not directly use the variables identified in the independent correlates analysis as training features. Instead, we used the Pearson correlation coefficient and a 10×10-fold cross-validated RFE process to identify the optimal feature set for each model from all available variables. During model development and validation, the 7 ML models exhibited distinct performance characteristics. In the training set, performance was evaluated using 10×10-fold internal cross-validation. Model discrimination was assessed using AUC, sensitivity, and specificity, whereas calibration was preliminarily evaluated using the Brier score. The K-Nearest Neighbors model demonstrated relatively poor discriminatory performance, with an AUC below the acceptable threshold of 0.7. Among the remaining models, XGBoost exhibited slightly superior discriminatory performance and marginally better calibration. In the test set, in addition to AUC, sensitivity, specificity, and Brier score, we evaluated model performance using calibration slope, calibration intercept, calibration curves, and DCA. Although the XGBoost model showed slightly better AUC, Brier score, calibration slope, and calibration intercept values than the RF model, these differences were not statistically significant. The RF model exhibited a slightly narrower CI for the Brier score and yielded the highest net benefit in DCA, indicating superior calibration stability and greater clinical utility. Based on findings from both the training and test sets, we believe this may be attributable to the RF model generating a smoother probability distribution, with sensitivity and specificity exhibiting less fluctuation across different decision thresholds. As the intended application of this model is to assist children’s guardians in achieving greater clinical benefit through risk assessment, overall clinical utility is particularly important. Accordingly, the RF model demonstrated stable performance and greater practical value, and was ultimately selected as the optimal model. ML models are often regarded as “black boxes” because of their limited transparency and explainability, which can reduce users’ acceptance of and trust in their predictions [33]. In recent years, the demand for Explainable Artificial Intelligence (XAI), a subfield of artificial intelligence, has grown substantially. XAI aims to improve the transparency, interpretability, and accountability of model behavior. One strength of this study is the use of SHAP for model interpretation, which has become one of the most widely used XAI methods in recent years. Compared with other interpretation methods, such as Local Interpretable Model-agnostic Explanations and Gradient-weighted Class Activation Mapping, SHAP offers the advantage of providing both global and local interpretability and is better suited for tree-based models such as RF and XGBoost [34]. In this study, SHAP was used to effectively illustrate the functionality of the RF model and explain how the RF model evaluates children based on features contained in the input data. The SHAP results indicated that children who are restless, speak loudly, have nonsoft muscles, and exhibit a sallow complexion are assigned a higher probability of current GHRS by the RF model. These findings are consistent with the independent correlates analysis and reflect characteristics that guardians can readily assess through repeated observation. Compared with similar studies, the explainability methods used in this research are relatively limited [34,35]. This approach was adopted because the primary users of the model and the intended audience for its explainability analyses are children’s guardians. Considering the nonmedical background of the target users, detailed interpretation results may be relatively difficult to understand and may offer greater reference value than practical utility. Therefore, rather than incorporating additional explainability methods, we included probability distribution plots to visually demonstrate the consistency between the disease probabilities predicted by the RF model and the observed trend in positive rates, thereby enhancing confidence in the model’s output. Furthermore, using the Streamlit framework, we integrated the RF model into a user-friendly online assessment platform. The platform collects the required evaluation information through a series of single-choice questions, with validation prompts displayed for unanswered items to ensure complete and structured data entry and to prevent logical inconsistencies. All questions are completed by the child’s guardian based on subjective judgment, consistent with the data collection approach used during model training, making the platform suitable for widespread use in routine health monitoring settings. Upon submission of the completed questionnaire, the platform displays the estimated probability of GHRS and the corresponding risk classification. Lifestyle recommendations are also provided to guide adjustments to the child’s daily habits and diet. These features demonstrate the feasibility of the model for large-scale dissemination and practical application.

The limitations of this study primarily include the following. First, to protect participant privacy and respect respondents’ preferences, all questionnaire items were optional, resulting in a substantial number of samples with missing values. As most features in this study were binary categorical variables, missing-value imputation produced clinically implausible contradictions among certain variables. To ensure data completeness and clinical relevance, we were required to exclude more than half of the samples containing missing values. Although the final sample size still substantially exceeded the required sample size, this approach resulted in considerable data loss and may have introduced some degree of bias into the findings. Second, although the model was developed using data from children in Longgang District, Shenzhen, China, and its performance was evaluated through internal cross-validation and validation in an independent test set, the test set samples were collected during the same period and from the same geographic region. In the absence of external validation across different periods or regions, there is insufficient evidence to confirm the model’s generalizability to other populations, despite the advantage of a large sample size. Third, all features used to develop the model consisted of sociodemographic information and data readily available in daily life. Although this enhances the model’s potential for widespread use in routine settings, the limited predictive relevance of these features suggests that there remains room for improvement in discriminatory performance. Fourth, the prevalence probability assessment models developed in this study were trained exclusively using traditional ML approaches. No reinforcement learning techniques were applied to optimize model performance, which may represent a methodological limitation [20]. Fifth, the use of SMOTETomek warrants cautious interpretation. As the outcome classes in this study were nearly balanced, resampling was not necessary to address severe class imbalance. Although sensitivity analyses demonstrated no meaningful deterioration in model performance, synthetic samples may still influence probability calibration. Future external validation studies and prospective evaluations should further compare the performance of simplified model pipelines that do not incorporate resampling. Finally, an important limitation is that GHRS status was determined using a self-rated diagnostic scale developed by the research team and completed by children’s guardians, rather than using an independent clinical reference standard. Although items from the diagnostic scale were excluded from model training, both the outcome and candidate predictors were derived from the same questionnaire source. This may have introduced shared-method bias and a degree of circularity in construct definition. Therefore, this model should be interpreted as estimating the probability that a child meets the study-specific GHRS scale criteria, rather than as an independent diagnostic tool for a clinically validated disease entity. To date, our team has made meaningful progress in GHRS-related research, including its definition, clinical manifestations, associated diseases, diagnostic approaches, pathogenic mechanisms, and risk factors, thereby establishing a preliminary research framework. Future efforts should focus on identifying and incorporating features more strongly associated with GHRS, expanding datasets to include samples from diverse geographic regions and periods, and integrating advanced modeling approaches, such as ensemble boosting techniques, into model development and validation. These improvements may further enhance the model’s predictive accuracy and clinical utility. Despite these limitations, the assessment model developed in this study remains practically useful for household self-assessment.

In summary, we successfully developed an explainable ML model to estimate the prevalence probability of GHRS in children using sociodemographic and clinical questionnaire data. During internal validation, the RF model demonstrated moderate discriminatory performance and acceptable calibration. It may serve as an exploratory tool for guardians to assess the likelihood that a child meets the study-specific GHRS scale criteria. However, external validation, usability testing, and prospective evaluation are required before implementation in clinical, community, or public health settings. In addition, this study differs from traditional approaches that either develop diagnostic models based on disease-specific symptoms using small cross-sectional samples or construct future disease risk prediction models using longitudinal cohort data. Instead, we used a large cross-sectional dataset to develop a GHRS prevalence probability assessment model based on sociodemographic characteristics and observable daily manifestations beyond disease-specific signs and symptoms. The predicted probability generated by the model reflects association-based estimates derived from the cross-sectional dataset. Models developed using this approach reduce reliance on professionals and specialized equipment for assessing disease-specific symptoms while quantifying disease probability using readily obtainable sociodemographic and behavioral information. As such, they may support exploratory self-assessment and hypothesis-generating population-level evaluations in settings involving nonprofessional users. Nevertheless, the practical screening value of this tool requires further confirmation through external validation, usability testing, and prospective studies. Future research should include external validation across different geographic regions, periods, and independent populations, as well as evaluation of the model’s usability and real-world impact in household and primary health care settings.